Coordinating on-demand service supply and demand paper notes

library(tidyverse)
library(triangle)
library(patchwork)
library(plotly)

plot_shaded_normal <- function(mean, sd, p = 0.5) {
  shade_up_to_x <- qnorm(p, mean = mean, sd = sd)

  data_normal <- data.frame(
    x = seq(mean - 4 * sd, mean + 4 * sd, length.out = 300)
  ) |>
    mutate(y = dnorm(x, mean = mean, sd = sd))

  shaded_data_normal <- data_normal |>
    filter(x <= shade_up_to_x)

  ggplot(data_normal, aes(x = x, y = y)) +
    geom_line() +
    geom_area(
      data = shaded_data_normal,
      aes(x = x, y = y),
      fill = "#ff8800",
      alpha = 0.3
    ) +
    geom_vline(xintercept = shade_up_to_x)
}

plot_shaded_triangle <- function(a, b, c = (a + b) / 2, p = 0.5) {
  shade_up_to_x <- triangle::qtriangle(p, a = a, b = b, c = c)

  data_tri <- data.frame(
    x = seq(a, b, length.out = 300)
  ) |>
    mutate(y = triangle::dtriangle(x, a = a, b = b, c = c))

  shaded_data_tri <- data_tri |>
    filter(x <= shade_up_to_x)

  ggplot(data_tri, aes(x = x, y = y)) +
    geom_line() +
    geom_area(
      data = shaded_data_tri,
      aes(x = x, y = y),
      fill = "#ff8800",
      alpha = 0.3
    ) +
    geom_vline(xintercept = shade_up_to_x)
}

plot_shaded_unif <- function(a = 0, b = 1, p = 0.5) {
  shade_up_to_x <- qunif(p, min = a, max = b)

  data_unif <- data.frame(
    x = seq(a, b, length.out = 300)
  ) |>
    mutate(y = dunif(x, min = a, max = b))

  shaded_data_unif <- data_unif |>
    filter(x <= shade_up_to_x)

  ggplot(data_unif, aes(x = x, y = y)) +
    geom_line() +
    geom_area(
      data = shaded_data_unif,
      aes(x = x, y = y),
      fill = "#ff8800",
      alpha = 0.3
    ) +
    geom_vline(xintercept = shade_up_to_x)
}

Reading paper:

Bai, J., So, K. C., Tang, C. S., Chen, X., & Wang, H. (2019). Coordinating supply and demand on an on-demand service platform with impatient customers. Manufacturing & Service Operations Management, 21(3), 556-570. https://doi.org/10.1287/msom.2018.0707

---

Slides

The slides used to be embedded… but they are too large to embed in this document anymore. A pdf version of the slides can be downloaded and viewed by clicking Here: Download Slides (PDF)

Notes

General overview and goals

• Customers arrive randomly
• Each service consists of random amount of service units (e.g. travel distance)
  • Assumed any service can be rendered by any agent (maybe not true in reality where drivers cancel trips based on distance)
• Fixed price rate \(p\) per service unit (e.g. dollars per kg)
• Fixed wage rate \(w\) per service unit
• Payout ratio: \(\frac{w}{p}\)

No focus on surge pricing for ease of modeling and cites customer dissatisfaction if price jumps in short time span.

Assumed \(p\) and \(w\) are known to customers and agents in advance.

Goal of service platform is to set \(p^*\) and \(w^*\) to maximize average profit.

Realized customer rate \(\lambda\) and price rate \(p\)

Setup

• Let \(\lambda\) be the customer request rate with max amount of service in a time period is \(\bar{\lambda}\)

• Customer value per service unit is \(v\) (heterogeneous with CDF \(F(\cdot)\))

• For a customer with valuation \(v\) and service request of \(D\) units the surplus is \((v - p)D\)

  • \(D\) is assumed to be independent of customer-type \(v\)

  • expected value of service units will be modeled as \(d = E(D)\)

• Utility function of a customer can be modeled as \(U(v) = (v - p)d - cW_q\)

  • \(v\) is value per service unit

  • \(p\) is price per service unit

  • \(d\) is expected value of service units

  • \(c\) is cost of waiting

  • \(W_q\) is expected wait time for service

• Realized customer request rate \(\lambda\) is the percent of \(\bar{\lambda}\) of customers who have \(U(v) \ge 0\) for a given \(p\)

  • aka \(\lambda = \text{Prob} \left\{ U(v) \ge 0 \right\} \cdot \bar{\lambda} \to \lambda = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\} \cdot \bar{\lambda}\)

  • steps

    • \((v-p)d-cW_q \ge 0\)

    • \(vd-pd-cW_q \ge 0\)

    • \(vd \ge pd-cW_q\)

    • \(v \ge p - \frac{c}{d}W_q\)

• This probability is labeled \(s\) (i.e. \(s = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\}\) and \(\lambda = s \bar{\lambda}\)) and represents the desired service rate. We can manipulate \(p\) to decide what percentage of \(\bar{\lambda}\) to target

• Since \(v \sim F(\cdot)\) then the price rate satisfies \(p = F^{-1} \left( 1 - \frac{\lambda}{\bar{\lambda}} \right) - \frac{c}{d} W_q\) (note this could be written \(p = F^{-1} \left( 1 - s \right) - \frac{c}{d} W_q\) to show service rate more directly)

  • Steps

    • \(\lambda = \text{Prob} \left\{ v \ge p + \frac{c}{d} W_q \right\} \cdot \bar{\lambda}\)

    • \(\lambda = \left[ 1 - F\left( p + \frac{c}{d} W_q \right) \right] \cdot \bar{\lambda}\)

    • \(\frac{\lambda}{\bar{\lambda}} = 1 - F\left( p + \frac{c}{d} W_q \right)\)

    • \(F\left( p + \frac{c}{d} W_q \right) = 1 - \frac{\lambda}{\bar{\lambda}}\)

    • \(p + \frac{c}{d} W_q = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right)\)

    • \(p = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\) or \(p = F^{-1}\left(1 - s\right) - \frac{c}{d} W_q\)

  • Holding everything else constant

    • Price rate decreases as wait time ( \(W_q\) ) increases

    • Price rate decreases as unit wait cost ( \(c\) ) increases

    • Price rate increases as number of service units ( \(d\) ) increases

Final results

Price is prescribed as:

\[ p = F^{-1}\left(1 - s\right) - \frac{c}{d} W_q = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \]

For given

• Service rate ( \(s = \frac{\lambda}{\bar{\lambda}}\) )
  • \(\lambda\) is requesting customers
  • \(\bar{\lambda}\) is max potential requesting customers
• Wait time ( \(W_q\) )
• Expected service units ( \(d\) )
• Cost of wait per unit time ( \(c\) )
• Inverse CDF of customer valuation of service per unit ( \(v \sim F(\cdot)\) )

# Price rate as a function of service level
price_rate <- function(lambda, lambda_bar, c, d, Wq, Finv = \(p) qunif(p = p, min = 0, max = 1)) {
  s <- lambda / lambda_bar
  Finv(1 - s) - (c / d) * Wq
}

How price varies with individual params

Holding all other factors constant and just manipulating the input shown on x axis.

These relationships are the similar in behavior regardless of choice of value distribution family ( \(F(\cdot)\) ).

• If \(c\) increases (and target service rate is the same), we charge less since the cost of waiting per unit of service \(\frac{c}{d}\) increases (if we kept same price we would be over priced and we’d be below target service)
• If \(d\) increases (and target service rate is the same), we charge more since the cost of waiting per unit of service \(\frac{c}{d}\) decreases (if we kept same price we would be under priced and we’d be above target service rate)
• If \(\lambda\) increases (aka target service rate increases), we charge less to mitigate the effects of waiting while maintaining value
• If \(W_q\) increases, we charge less to mitigate the effects of waiting while maintaining value

These notes don’t reflect how \(\lambda\) and \(W_q\) are related (i.e. more requests would lead to more wait time assuming the same supply).

ex_range <- 0:10

lambda_bar <- 10

lambda <- 2
c <- 2
d <- 2
Wq <- 2

# value dist
Finv <- \(p) qunif(p = p, min = 0, max = 1)


ex_lambda <- price_rate(ex_range, lambda_bar, ex_range, d, Wq, Finv)
ex_c <- price_rate(lambda, lambda_bar, ex_range, d, Wq, Finv)
ex_d <- price_rate(lambda, lambda_bar, c, ex_range, Wq, Finv)
ex_Wq <- price_rate(lambda, lambda_bar, c, d, ex_range, Finv)

func_behavior_df <- data.frame(
  ex_range = ex_range,
  lambda = ex_lambda,
  c = ex_c,
  d = ex_d,
  Wq = ex_Wq
)

func_behavior_df |>
  pivot_longer(-ex_range) |>
  mutate(name = ifelse(name == "c", "Cost of waiting (c)", name)) |>
  mutate(name = ifelse(name == "d", "Average service units (d)", name)) |>
  mutate(name = ifelse(name == "Wq", "Wait time (Wq)", name)) |>
  mutate(name = ifelse(name == "lambda", "Customer requests (lambda)", name)) |>
  ggplot(aes(x = ex_range, y = value)) +
  geom_line() +
  # geom_hline(yintercept = 0, linetype = "dotted") +
  facet_wrap(~name) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  theme(
    axis.text.y = element_blank(), axis.ticks.y = element_blank()
  ) +
  labs(
    y = "price", x = "",
    title = bquote("Behavior of price with inputs; " ~
      p == F^
        {
          -1
        } * (1 - lambda / bar(lambda) - frac(c, d) * W[q])),
    subtitle = bquote("Holding all other factors constant arbitrarily at 2 & " ~ bar(lambda) == 10)
  )

Uniform distribution of value \(v \sim F_{unif_{[0, 1]}}(\cdot)\)

The paper will use a uniform distribution \([0, 1]\) to simplify analysis but insights extend to other families

c <- 0.25
d <- 2
Wq <- 0.5

# value dist (unifo) params
Finv <- \(p) qunif(p = p, min = 0, max = 1)

lambda <- 2
lambda_bar <- 10
s <- lambda / lambda_bar
p1 <- plot_shaded_unif(a = 0, b = 1, p = s) +
  labs(
    title = sprintf(
      "For service level %.0f%% -> price = %.2f ",
      s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
    ),
    subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )


lambda <- 6
lambda_bar <- 10
s <- lambda / lambda_bar
p2 <- plot_shaded_unif(a = 0, b = 1, p = s) +
  labs(
    title = sprintf(
      "For service level %.0f%% -> price = %.2f ",
      s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
    ),
    subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )


p1 + p2

Normal distribution of value \(v \sim F_{normal}(\cdot)\)

c <- 2
d <- 2
Wq <- 2

# value dist (normal) params
mu <- 5
sd <- 1
Finv <- \(p) qnorm(p = p, mean = mu, sd = sd)

lambda <- 3
lambda_bar <- 10
s <- lambda / lambda_bar
p1 <- plot_shaded_normal(mean = mu, sd = sd, p = s) +
  labs(
    title = sprintf(
      "For service level %.0f%% -> price = %.2f ",
      s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
    ),
    subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

lambda <- 8
lambda_bar <- 10
s <- lambda / lambda_bar
p2 <- plot_shaded_normal(mean = mu, sd = sd, p = s) +
  labs(
    title = sprintf(
      "For service level %.0f%% -> price = %.2f ",
      s * 100, price_rate(lambda, lambda_bar, c, d, Wq, Finv)
    ),
    subtitle = sprintf("Given c = %.1f; d = %.1f, Wq = %.1f", c, d, Wq),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

p1 + p2

Realized providers \(k\) and wage rate \(w\)

Setup

• \(K\) - maximum number of potential providers
• \(k\) - realized number of providers for a given price \(p\) and wage rate \(w\) (wage as a percentage - \(\alpha\) - of \(p\))
  • \(k \le K\)
• \(\mu\) - average service speed; so \(\frac{\mu}{d}\) is service rate (i.e. average number of customers served per hour; \(d\) is average service units requested)
• Given realized customers ( \(\lambda\) ) and realized providers ( \(k\) ), the utilization is \(\frac{\lambda}{k \cdot (\mu / d)} = \frac{\lambda d}{k \mu}\)
  • Aka \(\rho = \frac{\lambda d}{k \mu}\) for \(M/M/k\) queue
  • Similar to typical queueing \(\rho = \frac{\lambda}{\mu}\), but incorporating \(k\) providers and \(d\) average service units
• Wage per unit time of a single provider is \(w\mu\)
• Average earning rate is \(w\frac{\lambda d}{k}\)
  • that is: single provider earning \(w \mu\) times utilization \(\frac{\lambda d}{k \mu}\)
  • \(w \mu \cdot \frac{\lambda d}{k \mu} = w\frac{\lambda d}{k}\)
• Each provider has a “reservation earning rate” \(r\) that has CDF \(G(\cdot)\)
  • i.e. they won’t average earning rate \(\ge\) than reservation rate (aka \(w\frac{\lambda d}{k} \ge r\) )
• \(\beta\) will denote proportion of providers partipating
  • \(\beta = \text{Prob} \left\{ r \le w(\lambda d / k) \right\} = G \left( w(\lambda d / k) \right)\)
  • Also: \(G^{-1}(\beta) = w \frac{\lambda d}{k}\)
  • This means realized providers \(k\) is \(k = \beta K\)

Final results

Wage is prescribed as:

\[ w = G^{-1}(\beta) \frac{k}{\lambda d} = G^{-1}(\frac{k}{K}) \frac{k}{\lambda d} \]

For given

• Participation rate ( \(\beta = \frac{k}{K}\) )
  • \(k\) is participating providers
  • \(K\) is max potential participating providers
• Expected service units ( \(d\) )
• Number of requesting customers ( \(\lambda\) )
• Inverse CDF of provider reservation rate ( \(r \sim G(\cdot)\) )

# Wage rate as a function of provider participation rate
wage_rate <- function(k, K, lambda, d, Ginv = \(p) qunif(p = p, min = 0, max = 1)) {
  Ginv(k / K) * (k / (lambda * d))
}

How wage varies with individual params

Holding all other factors constant and just manipulating the input shown on x axis.

These relationships are the similar in behavior regardless of choice of reservation rate distribution family ( \(G(\cdot)\) ).

• If \(\lambda\) increases (and target participation rate is the same), we pay less since the average utilization of providers increases (more fares to go around)
• If \(d\) increases (and target participation rate is the same), we pay less since the pay out per fare increases (more units of pay)
• If \(k\) increases (aka target participation rate increases), we need to pay more to incentivize providers to enter

ex_range <- 0:10

K <- 10

k <- 2
lambda <- 2
d <- 2

# reservation rate dist
Ginv <- \(p) qunif(p = p, min = 0, max = 1)

ex_k <- wage_rate(ex_range, K, lambda, d, Ginv)
ex_lambda <- wage_rate(k, K, ex_range, d, Ginv)
ex_d <- wage_rate(k, K, lambda, ex_range, Ginv)

func_behavior_df <- data.frame(
  ex_range = ex_range,
  k = ex_k,
  lambda = ex_lambda,
  d = ex_d
)

func_behavior_df |>
  pivot_longer(-ex_range) |>
  mutate(name = ifelse(name == "k", "Pariticpating providers (k)", name)) |>
  mutate(name = ifelse(name == "d", "Average service units (d)", name)) |>
  mutate(name = ifelse(name == "lambda", "Customer requests (lambda)", name)) |>
  ggplot(aes(x = ex_range, y = value)) +
  geom_line() +
  # geom_hline(yintercept = 0, linetype = "dotted") +
  facet_wrap(~name) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  theme(
    axis.text.y = element_blank(), axis.ticks.y = element_blank()
  ) +
  labs(
    y = "wage", x = "",
    title = bquote("Behavior of wage with inputs; " ~
      w == G^
        {
          -1
        } * (frac(k, K)) * frac(k, lambda * d)),
    subtitle = bquote("Holding all other factors constant arbitrarily at 2 & " ~ K == 10)
  )

Uniform distribution of value \(r \sim G_{unif_{[0, 1]}}(\cdot)\)

The paper will use a uniform distribution \([0, 1]\) to simplify analysis but insights extend to other families

lambda <- 2
d <- 2

# value dist (unifo) params
Ginv <- \(p) qunif(p = p, min = 0, max = 1)

k <- 2
K <- 10
beta <- k / K
p1 <- plot_shaded_unif(a = 0, b = 1, p = beta) +
  labs(
    title = sprintf(
      "For partipication level %.0f%% -> wage = %.2f ",
      beta * 100, wage_rate(k, K, lambda, d, Ginv)
    ),
    subtitle = sprintf("Given lambda = %.1f; d = %.1f", lambda, d),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )


k <- 8
K <- 10
beta <- k / K
p2 <- plot_shaded_unif(a = 0, b = 1, p = beta) +
  labs(
    title = sprintf(
      "For service level %.0f%% -> wage = %.2f ",
      beta * 100, wage_rate(k, K, lambda, d, Ginv)
    ),
    subtitle = sprintf("Given lambda = %.1f; d = %.1f", lambda, d),
    x = "v",
    y = ""
  ) +
  theme(
    plot.title = element_text(size = 10),
    plot.subtitle = element_text(size = 8),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )


p1 + p2

Objective: max profit

• Profit ( \(\pi\) ) is \(\pi = \lambda (p - w) d\)
  • \(\lambda\) - customer requests
  • \(p\) - price rate per service unit
  • \(w\) - wage rate per service unit
  • \(d\) - average service units
• We can plug in the definitions of \(p\) (in terms of service rate - \(\lambda\)) and \(w\) (in terms of participating provider rate - \(k\))
  • \(\pi(k, \lambda) = \lambda d \left[ \left( F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \right) - \left( G^{-1}(\frac{k}{K}) \frac{k}{\lambda d} \right) \right]\)
  • Or written with \(s\) to denote request rate and \(\beta\) to denote provider participation rate: \(\max_{k, \lambda} \pi(k, \lambda) \equiv \lambda d \left[ \left( F^{-1}\left(1 - s\right) - \frac{c}{d} W_q \right) - \left( G^{-1}(\beta) \frac{k}{\lambda d} \right) \right]\)
  • Subject to \(\frac{\lambda d}{k \mu} < 1\) (i.e. the requests don’t outpace the service)
    • Note: maybe not a realistic constraint for peak event times, but maybe surge pricing and wait times can diminish request rate enough to hold
• We’ll use \(\pi(k, \lambda)\) to find optimal supply \(k^*\) and optimal demand \(\lambda^*\), and in turn we can use \(k^*\) and \(\lambda^*\) to find optimal price \(p^*\) and optimal wage rate \(w^*\)

Base model

• Wage rate ( \(w\) ) will be a proportion ( \(\alpha\) ) of price rate ( \(p\) ) - \(w = \alpha p\)

• Waiting time \(W_q\) will be based on \(M/M/k\) queue with arrival rate \(\lambda\) and service rate \(\frac{\mu}{d}\)

  • \(W_q = \frac{1}{1 + \left( \frac{k! (1 - \rho)}{k^k \rho^k} \right) \sum_{i=0}^{k-1} \frac{k^i \rho^i}{i!}} \left[ \frac{\rho}{\lambda (1 - \rho)} \right]\)
  • Where \(\rho = \frac{\lambda d}{k \mu}\) is system utilization with \(\rho < 1\)

• For simplicity in base model:

  • The distribution of \(v\) (so far referred to as \(F(\cdot)\) ) and \(r\) (so far referred to as \(G(\cdot)\) ) will both be modeled as uniform on range \([0, 1]\). This means we can drop out mention of them from current model.

    • \(p = F^{-1}\left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\) becomes \(p = \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q\)

    • \(w = G^{-1}(\frac{k}{K}) \frac{k}{\lambda d}\) becomes \(w = \frac{k}{K} \cdot \frac{k}{\lambda d} = \frac{k^2}{K \lambda d}\)

    • \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left[ \left( F^{-1}\left(1 - s\right) - \frac{c}{d} W_q \right) - \left( G^{-1}(\beta) \frac{k}{\lambda d} \right) \right]\) becomes \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left[ \left( \left(1 - s \right) - \frac{c}{d} W_q \right) - \left( \frac{k^2}{K \lambda d} \right) \right]\)

  • Wage rate is \(\alpha\) proportion of \(p\) : \(w = \alpha p \to \frac{w}{\alpha} = p\) ; avoid direct \(p\) formulation since it becomes way messier when incorporating Wait time formulation

    • \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d (p - w) = \lambda d \left( \frac{w}{\alpha} - w \right)\)

    • Incorporating definition of \(w = \frac{k^2}{K \lambda d}\) we can write as \(\max_{k, \lambda} \pi(k, \lambda) = \lambda d \left( \frac{k^2}{K \lambda d \alpha} - \frac{k^2}{K \lambda d} \right) = \frac{k^2}{K \alpha} - \frac{k^2}{K} = \frac{k^2 - k^2 \alpha}{K \alpha} = \frac{k^2 (1 - \alpha)}{K \alpha}\)

    • Pretty clear to see we want to maximize number of participating providers \(k\)

  • This is subject to constraints

    • \(\frac{\lambda d}{k \mu} < 1\)

    • \(w = \alpha p\) which can be written as (still avoiding writing out \(W_q\) here; see def above)

      • \(\frac{k^2}{K \lambda d} = \alpha \left[ \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \right]\)

      • This implies \(k^2 = {K \lambda d}\alpha \left[ \left(1 - \frac{\lambda}{\bar{\lambda}}\right) - \frac{c}{d} W_q \right]\)

exact_wait_time_ <- function(k, lambda, mu, d) {
  # $$
  # W_q = \frac{1}{1 + \left( \frac{k! (1 - \rho)}{k^k \rho^k} \right) \sum_{i=0}^{k-1} \frac{k^i \rho^i}{i!}} \left[ \frac{\rho}{\lambda (1 - \rho)} \right]
  # $$

  # Calculate rho (system utilization)
  rho <- (lambda * d) / (k * mu)
  stopifnot(rho < 1)

  # First fraction term
  # ... Denominator term 1
  ## scalar to adjust probs of term 2
  term1 <- (factorial(k) * (1 - rho)) / (k^k * rho^k)

  # ... Denominator term 2
  ## each term in summation is prob of i servers busy
  i <- 0:(k - 1)
  term2 <- sum((k^i * rho^i) / factorial(i))

  # Full formula for Wq
  wait_time <- 1 / (1 + term1 * term2) * (rho / (lambda * (1 - rho)))
  return(wait_time)
}

exact_wait_time <- Vectorize(exact_wait_time_, c("k", "lambda", "mu", "d"))

prob_queue_len <- function(i, k, lambda, mu, d) {
  rho <- (lambda * d) / (k * mu)

  n <- 0:(k - 1)
  p_0 <- 1 / (sum(rho^n / factorial(n)) + (rho^k / factorial(k)) * (1 / (1 - rho)))

  if (i == 0) {
    p_i <- p_0
  } else if (i < k) {
    p_i <- (rho^i / factorial(i)) * p_0
  } else {
    p_i <- (rho^i / (factorial(k) * k^(i - k))) * p_0
  }

  return(p_i)
}

queue_len_probs <- function(k, lambda, mu, d) {
  probs <- c()
  i_vals <- 0:(k + 10)
  for (i in i_vals) {
    p_i <- prob_queue_len(i, k, lambda, mu, d)
    probs[i + 1] <- p_i
  }

  return(data.frame(len = i_vals, prob = probs))
}

check_rho_constraint <- function(lambda, d, k, mu) {
  rho <- (lambda * d) / (k * mu)

  return(rho < 1)
}

check_fixed_payout_constraint <- function(
    k, K, lambda, lambda_bar,
    c, d,
    mu,
    alpha,
    tolerance = 0.01) {
  v1 <- k^2

  Wq <- exact_wait_time(k, lambda, mu, d)
  v2 <- K * lambda * d * alpha * ((1 - (lambda / lambda_bar)) - (c / d) * Wq)

  return(isTRUE(all.equal(v1, v2, tolerance = tolerance)))
}

profit_function <- function(
    k, K, lambda, lambda_bar,
    c, d,
    mu,
    alpha,
    Finv = \(p) qunif(p, min = 0, max = 1),
    Ginv = \(p) qunif(p, min = 0, max = 1)) {
  meets_rho_constraint <- check_rho_constraint(lambda, d, k, mu)

  if (!meets_rho_constraint) {
    return(list(
      meets_rho_constraint = FALSE,
      meets_fixed_payout_constraint = FALSE
    ))
  }

  meets_fixed_payout_constraint <- check_fixed_payout_constraint(k, K, lambda, lambda_bar, c, d, mu, alpha)

  if (!meets_fixed_payout_constraint) {
    return(list(
      meets_rho_constraint = TRUE,
      meets_fixed_payout_constraint = FALSE
    ))
  }

  Wq <- exact_wait_time(k, lambda, mu, d)

  price <- price_rate(lambda, lambda_bar, c, d, Wq, Finv)
  wage <- wage_rate(k, K, lambda, d, Ginv)

  profit <- (k^2 * (1 - alpha)) / (K * alpha)

  list(
    k = k,
    K = K,
    lambda = lambda,
    lambda_bar = lambda_bar,
    c = c,
    d = d,
    mu = mu,
    rho = (lambda * d) / (k * mu),
    alpha = alpha,
    Wq = Wq,
    Lq = lambda * Wq,
    queue_len_dist = queue_len_probs(k, lambda, mu, d),
    price = price,
    wage = wage,
    profit = profit,
    meets_rho_constraint = meets_rho_constraint,
    meets_fixed_payout_constraint = meets_fixed_payout_constraint
  )
}

display_results <- function(res) {
  display_res <- res[c("k", "K", "lambda", "lambda_bar", "rho", "Wq", "price", "wage", "alpha", "profit")]

  display_res |>
    as.data.frame() |>
    t() |>
    as.data.frame() |>
    rownames_to_column() |>
    set_names(c("Param", "Value")) |>
    mutate(Value = as.numeric(sprintf("%.3f", Value)))
}

Here is the results for row 1 of Table 1 to show my formulation matches.

res <- profit_function(k = 7, K = 50, lambda = 2.71, lambda_bar = 10, c = 1, d = 1, mu = 1, alpha = 0.5)

res$queue_len_dist |>
  filter(prob > 0.01) |>
  ggplot(aes(x = len, y = prob)) +
  geom_bar(stat = "identity")

display_results(res)

Interestingly, and unreported by the paper. This is not a unique solution. Their Table 1 seems to prefer a solution with higher price and lower wait times.

res <- profit_function(k = 7, K = 50, lambda = 4.73, lambda_bar = 10, c = 1, d = 1, mu = 1, alpha = 0.5)

res$queue_len_dist |>
  filter(prob > 0.01) |>
  ggplot(aes(x = len, y = prob)) +
  geom_bar(stat = "identity")

display_results(res)

Examples searching for max profit

My formulation matches paper results, but to save time I’ve used a lower precision. Because of this, my optimal values for lambda will be off from paper table by decimal amounts.

find_optimal_k_and_lambda <- function(
    K, lambda_bar, c, d, mu, alpha,
    Finv = \(p) qunif(p, min = 0, max = 1),
    Ginv = \(p) qunif(p, min = 0, max = 1),
    prefer_smaller_price = TRUE,
    k_step = 1,
    lambda_step = 0.01,
    tied_profit_tolerance = 0.001) {
  best_outcome <- NULL

  possible_k <- seq(1, K, by = k_step)

  # if (interactive()) {
  #   n_iter <- length(possible_k) * length(possible_lambda)
  #   pb <- txtProgressBar(min = 0, max = n_iter, style = 3)
  #   progress <- 0
  # }

  for (k in possible_k) {
    possible_lambda <- seq(lambda_step, k, by = lambda_step)
    for (lambda in possible_lambda) {
      res <- profit_function(
        k = k, K = K,
        lambda = lambda, lambda_bar = lambda_bar,
        c = c, d = d,
        mu = mu,
        alpha = alpha,
        Finv = Finv, Ginv = Ginv
      )

      if (res$meets_rho_constraint & res$meets_fixed_payout_constraint) {
        if (is.null(best_outcome)) {
          best_outcome <- res
          next
        }

        if (res$profit > best_outcome$profit) {
          best_outcome <- res
          next
        }

        # if profit is tied with best (within a tolerance)
        if (isTRUE(all.equal(res$profit, best_outcome$profit, tolerance = tied_profit_tolerance))) {
          if (prefer_smaller_price & res$price < best_outcome$price) {
            best_outcome <- res
          } else if (!prefer_smaller_price & res$price > best_outcome$price) {
            best_outcome <- res
          }
        }
      }

      # if (interactive()) {
      #   progress <- progress + 1
      #   setTxtProgressBar(pb, progress)
      # }
    }
  }

  # if (interactive()) {
  #   close(pb)
  # }
  return(best_outcome)
}

Here is the results for row 1 of Table 1 to show my formulation can find the same solution, but with a lack of precision.

res <- find_optimal_k_and_lambda(
  K = 50,
  lambda_bar = 10,
  c = 1,
  d = 1,
  mu = 1,
  alpha = 0.5,
  prefer_smaller_price = FALSE
)

display_results(res)

res <- find_optimal_k_and_lambda(
  K = 50,
  lambda_bar = 10,
  c = 1,
  d = 1,
  mu = 1,
  alpha = 0.5,
  prefer_smaller_price = TRUE
)

display_results(res)

How inputs affect shape (integer k)

\(\alpha\) - wage rate percentage of price

Preferring higher prices

all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = 10,
    c = 1,
    d = 1,
    mu = 1,
    alpha = alpha,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
  pivot_longer(-alpha) |>
  ggplot(aes(x = alpha, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

Preferring lower prices

all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = 10,
    c = 1,
    d = 1,
    mu = 1,
    alpha = alpha,
    prefer_smaller_price = TRUE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
  pivot_longer(-alpha) |>
  ggplot(aes(x = alpha, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

\(\bar{\lambda}\) - max requests

Preferring higher prices

all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = lambda_bar,
    c = 1,
    d = 1,
    mu = 1,
    alpha = 0.5,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
  pivot_longer(-lambda_bar) |>
  ggplot(aes(x = lambda_bar, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

Preferring lower prices

all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = lambda_bar,
    c = 1,
    d = 1,
    mu = 1,
    alpha = 0.5,
    prefer_smaller_price = TRUE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
  pivot_longer(-lambda_bar) |>
  ggplot(aes(x = lambda_bar, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

\(d\) - average service units

Preferring higher prices

all_res <- list()
i <- 0
for (d in seq(1, 10, 2)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = 10,
    c = 1,
    d = d,
    mu = 1,
    alpha = 0.5,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, d, rho, Wq, price, wage, profit) |>
  pivot_longer(-d) |>
  ggplot(aes(x = d, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

Preferring lower prices

all_res <- list()
i <- 0
for (d in seq(1, 10, 2)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = 10,
    c = 1,
    d = d,
    mu = 1,
    alpha = 0.5,
    prefer_smaller_price = TRUE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, d, rho, Wq, price, wage, profit) |>
  pivot_longer(-d) |>
  ggplot(aes(x = d, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

How inputs affect shape (non-integer k)

\(\alpha\)

all_res <- list()
i <- 0
for (alpha in seq(0.1, 0.9, 0.1)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = 20,
    c = 1,
    d = 1,
    mu = 1,
    alpha = alpha,
    k_step = 0.2,
    lambda_step = 0.2,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, rho, alpha, Wq, price, wage, profit) |>
  pivot_longer(-alpha) |>
  ggplot(aes(x = alpha, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

\(\bar{\lambda}\)

all_res <- list()
i <- 0
for (lambda_bar in seq(10, 100, 10)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = 50,
    lambda_bar = lambda_bar,
    c = 1,
    d = 1,
    mu = 1,
    alpha = 0.5,
    k_step = 0.2,
    lambda_step = 0.2,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, lambda_bar, rho, Wq, price, wage, profit) |>
  pivot_longer(-lambda_bar) |>
  ggplot(aes(x = lambda_bar, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

\(K\)

all_res <- list()
i <- 0
for (K in seq(10, 100, 10)) {
  i <- i + 1

  res <- find_optimal_k_and_lambda(
    K = K,
    lambda_bar = 20,
    c = 1,
    d = 1,
    mu = 1,
    alpha = 0.5,
    k_step = 0.2,
    lambda_step = 0.2,
    prefer_smaller_price = FALSE
  )
  res$queue_len_dist <- NULL
  all_res[[i]] <- res
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

res_df |>
  select(k, lambda, K, rho, Wq, price, wage, profit) |>
  pivot_longer(-K) |>
  ggplot(aes(x = K, y = value)) +
  geom_line() +
  geom_point() +
  facet_wrap(~name, scales = "free_y")

Varying \(K\) and \(\bar{\lambda}\)

all_res <- list()
i <- 0
for (K in seq(10, 100, 20)) {
  for (lambda_bar in seq(10, 100, 20)) {
    i <- i + 1

    res <- find_optimal_k_and_lambda(
      K = K,
      lambda_bar = lambda_bar,
      c = 1,
      d = 1,
      mu = 1,
      alpha = 0.5,
      k_step = 1,
      lambda_step = 1,
      prefer_smaller_price = FALSE
    )
    res$queue_len_dist <- NULL
    all_res[[i]] <- res
  }
}

res_df <- do.call(rbind, lapply(all_res, as.data.frame))

res_df

plot_cols <- c("k", "lambda", "rho", "Wq", "price", "wage", "profit")
res_df |>
  mutate_at(all_of(plot_cols), scale) |>
  select(K, lambda_bar, k, lambda, rho, Wq, price, wage, profit) |>
  pivot_longer(all_of(plot_cols)) |>
  rename(`z-score` = value) |>
  ggplot(aes(x = K, y = lambda_bar, fill = `z-score`)) +
  geom_tile() +
  facet_wrap(~name)

An approximation scheme

• To approximate and give tractable analytical results can instead be approximated with \(W_q = \frac{ \rho ^ \sqrt{2 (k + 1)} }{ \lambda (1 - \rho) }\)

  • The approximation formula is exact for an \(M/M/1\) queue (i.e. plug \(k = 1\) into above complex form and it reduces to \(W_q = \frac{ \rho ^ 2 }{ \lambda (1 - \rho) } = \frac{ \rho ^ \sqrt{2 (k + 1)} }{ \lambda (1 - \rho) }\) )
  • The approximation has been shown (see Sakasegawa 1977) to provide a very good estimate when \(k > 1\) .

See below for how this approximation compares to more exact, but more unwieldy formulation.

The paper actually goes further from exact calc and rewrites to be \(W_q = \frac{ \rho ^ \sqrt{2 (n + 1)} }{ \lambda (1 - \rho) }\). And \(n\) is set to be \(k^*\). This means for \(k^*\), \(W_q\) is well approximated, but for values of \(k < k^*\) the \(W_q\) will be over estimated and for values of \(k > k^*\), \(W_q\) will be underestimated.

approx_wait_time_ <- function(k, lambda, mu, d) {
  # $$
  # W_q = \frac{ \rho ^ \sqrt{2 (k + 1)} }{ \lambda (1 - \rho) }
  # $$

  # Calculate rho (system utilization)
  rho <- (lambda * d) / (k * mu)

  # numerator shows traffic intensity - for large k, less sensitive to traffic
  # denominator combines arrival rate and idle capacity
  # to show traffic that can be handled without queueing
  wait_time <- (rho^sqrt(2 * (k + 1))) / (lambda * (1 - rho))
  return(wait_time)
}

approx_wait_time <- Vectorize(approx_wait_time_, c("k", "lambda", "mu", "d"))

k <- 1:50

lambda <- 4
mu <- 16
d <- 2

wait_time_compare <- data.frame(
  k = k,
  Exact = exact_wait_time(k, lambda, mu, d),
  Approx = approx_wait_time(k, lambda, mu, d)
) |>
  pivot_longer(-k, names_to = "Formula") |>
  mutate(`log(Wait time)` = log(value)) |>
  rename(`Wait time` = value) |>
  pivot_longer(contains("Wait"), names_to = "trans") |>
  mutate(trans = forcats::fct_rev(trans), Formula = forcats::fct_rev(Formula))

ggplot(wait_time_compare, aes(x = k, y = value, color = Formula, linetype = Formula)) +
  geom_line() +
  facet_wrap(~trans, scales = "free_y") +
  labs(
    title = "Comparison of Approx and Exact M/M/k wait time formulas",
    subtitle = bquote(lambda == .(lambda) ~ ", " ~
      mu == .(mu) ~ ", " ~
      d == .(d) ~ " --> " ~
      rho == frac(lambda * d, k * mu) ~ ">=" ~ 1)
  )

Finding optimal alpha for given \(K\) & \(\bar{\lambda}\)

This is a remake of Table 6 in the paper to apply conditional formatting.

# Example data: Replace this with your actual data
data <- data.frame(
  lambda = seq(10, 100, by = 10), # Row values (lambda)
  K10 = c(0.67, 0.67, 0.67, 0.67, 0.67, 0.68, 0.68, 0.68, 0.68, 0.68),
  K20 = c(0.49, 0.52, 0.54, 0.55, 0.56, 0.56, 0.57, 0.57, 0.57, 0.57),
  K30 = c(0.42, 0.46, 0.48, 0.49, 0.51, 0.51, 0.52, 0.53, 0.53, 0.53),
  K40 = c(0.36, 0.41, 0.44, 0.46, 0.47, 0.48, 0.49, 0.49, 0.50, 0.50),
  K50 = c(0.33, 0.37, 0.40, 0.42, 0.44, 0.45, 0.46, 0.47, 0.48, 0.48),
  K60 = c(0.30, 0.34, 0.38, 0.40, 0.42, 0.43, 0.44, 0.45, 0.46, 0.46),
  K70 = c(0.27, 0.32, 0.35, 0.38, 0.39, 0.41, 0.42, 0.43, 0.44, 0.45),
  K80 = c(0.25, 0.30, 0.33, 0.36, 0.38, 0.39, 0.40, 0.41, 0.42, 0.43),
  K90 = c(0.24, 0.28, 0.31, 0.34, 0.36, 0.38, 0.39, 0.40, 0.41, 0.42),
  K100 = c(0.22, 0.27, 0.30, 0.32, 0.34, 0.36, 0.37, 0.39, 0.40, 0.41)
)

# Reshape data to long format
long_data <- data %>%
  pivot_longer(
    cols = starts_with("K"), # Select columns K10, K20, ..., K100
    names_to = "K",
    names_prefix = "K",
    values_to = "value"
  ) |>
  mutate(K = as.numeric(K))

# Create a heatmap
ggplot(long_data, aes(x = K, y = lambda, fill = value)) +
  geom_tile() +
  geom_text(aes(label = sprintf("%.2f", value)), color = "black", size = 3) +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(x = "K", y = expression(bar(lambda)), fill = expression(alpha)) +
  theme_minimal() +
  theme(
    panel.grid = element_blank(),
    axis.title.y = element_text(angle = 0, vjust = 0.5)
  ) +
  scale_y_reverse()